Ab initio molecular dynamics study of the desorption of D2 from Si (100) 
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Ab initio molecular dynamics calculations of deuterium desorbing from Si (100) have been performed in 
order to monitor the energy redistribution among the various D2 and silicon degrees of freedom during 
the desorption process. The calculations show that a considerable part of the potential energy at the 
transition state to desorption is transferred to the silicon lattice. The deuterium molecules leave the 
surface vibrationally hot and rotationally cold, in agreement with thermal desorption experiments; the 
mean kinetic energy, however, is larger than found in a laser-induced desorption experiment. We discuss 
possible reasons for this discrepancy. 
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Hydrogen adsorption on and desorption from Si sur- 
faces are of great technological relevance for, e.g., the 
etching and passivation of Si surfaces or the growth of Si 
crystals (see, e.g., Ref. [Q and references therein). Be- 
sides, the dynamics of the hydrogen interaction with Si 
surfaces is also of fundamental interest caused, among 
others, by the so-called barrier puzzle: While the stick- 
ing coefficient of molecular hydrogen on Si surfaces is 
very small ^|||],[| indicating a high barrier to adsorp- 
tion, in desorption experiments an almost thermal mean 
kinetic energy of the molecules was found Q indicating 
a low adsorption barrier. In order to explain this puzzle 
it was suggested to take the strong surface rearrange- 
ment of Si upon hydrogen adsorption into account ||: 
The hydrogen molecules impinging on the Si substrate 
from the gas phase typically encounter a Si configuration 
which is unfavorable for dissociation, while desorbing hy- 
drogen molecules leave the surface from a rearranged Si 
configuration with a low barrier. Unfortunately, at this 
point it is unclear how low the barrier really is, because 
the results for the dissociative adsorption probability of 
references [9 and || differ by almost three orders of mag- 
nitude. This difference translates into a difference of the 
minimum barrier heights measured in these two experi- 
ments of about 0.4 eV. 

Total-energy calculations using the cluster approach 
jj],||]!,|l(J have found activation barriers for associative 
desorption of about 1 eV. With the assumption of a lat- 
tice rearrangement energy of about 0.7 eV the experi- 
mental adsorption results of Refs. Q|J and the desorp- 
tion results of Ref. |J could be reproduced by quantum 
dynamical model calculations (l^Jl^]. Still the exact ad- 
sorption and desorption mechanism is strongly debated 
in the system H^/Si (100). In view of density-functional 
theory (DFT) calculations for Si (100) a lattice-relaxation 
energy of 0.7 eV seemed to be too high p3| . In- 
deed, in detailed DFT calculations of the H 2 /Si(100) 
potential energy surface (PES) using the supercell ap- 
proach Q,|l5],|l6| the adsorption barriers were found to 
be only 0.3 - 0.4 eV with a substrate rearrangement en- 
ergy of about 0.15 eV. 

Based on the slab calculations |L4j and approximating 
the high-dimensional PES by a three-dimensional one, 



two different quantum dynamical studies were performed 
|]l7| , ^8| , where the relaxation of the Si substrate upon hy- 
drogen adsorption was represented by one idealized Si 
phonon coordinate. Although both calculations used the 
same ab initio energies as a source, the results of the dy- 
namical calculations did not agree quantitatively. How- 
ever, in both studies the dynamical coupling of the des- 
orption path to the Si substrate was low. 

In order to help clarifying these various confusing and 
in fact partially conflicting results we have studied the 
associative desorption of D2 from Si (100) by determin- 
ing how the potential energy at the barrier is distributed 
over the various degrees of freedom of this system (hy- 
drogen vibration, rotation, and translational energy, and 
vibrations of the Si substrate). The calculation of PESs 
is an important prerequisite for understanding reaction 
dynamics. For a quantitative analysis, however, a calcu- 
lation of the dynamics is indispensable. Since the sub- 
strate relaxation plays such an important role for the 
hydrogen desorption from Si (100), a proper description 
of the Si dynamics during the desorption event is re- 
quired. We feel that a reliable dynamical calculation 
requires the consideration of the six degrees of freedom 
of the hydrogen molecule plus at least two degrees of 
freedom of the Si surface (dimer stretch and tilt). A 
quantum dynamical treatment of the desorption dynam- 
ics taking these eight degrees of freedom simultaneously 
into account is presently not possible. The record still 
stands at six dimensions jl9|]. We have therefore per- 
formed ab initio molecular dynamics calculations to mon- 
itor the energy distribution of D2 molecules desorbing 
from Si (100). This allows us to assess the dynamical 
consequences of the calculated PES. We will show that 
the deuterium molecules leave the surface vibrationally 
hot and rotationally cold, in agreement with thermal des- 
orption experiments pfj| ]. Indeed, approximately 0.1 eV 
of the potential energy at the transition state is trans- 
ferred to substrate vibrations. Still, also the molecules 
receive a noticeable kinetic energy which is at variance 
with a laser-induced desorption experiment 0. We will 
discuss possible reasons for this interesting discrepancy. 



In the ab initio molecular dynamics approach [ET 22 



the forces necessary to integrate the classical equations 
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FIG. 1. Snapshots of a trajectory of D2 desorbing from 
Si (100) starting at the transition state with the Si atoms ini- 
tially at rest. The dark Si atoms correspond to the Si positions 
after the desorption event. 



of motion are determined by DFT calculations. The 
exchange-correlation functional is treated in the general- 
ized gradient approximation (GGA) . For the hydro- 
gen atoms the full 1 jr potential is used. In previous slab 
studies the total energies were calculated within the local 
density approximation (LDA) with a posteriori GGA cor- 
rections The main effects of using the GGA in the 
complete self-consistent cycle are a small increase of the 
theoretical lattice constant of Si |03 and a slight rise in 
the barrier height from E b = 0.3 cV p| to E b = 0.4 eV. 
To correctly represent the up and down buckling of the 
clean Si(100) surface we use a (2x2) surface unit cell. 
The Si slab consists of five atomic layers. The topmost 
three of them are free to move in the molecular dynam- 
ics simulations, while the remaining two layers are fixed 
at their bulk positions. The density-functional calcula- 
tions are performed with two k-points in the irreducible 
part of the Brillouin zone and 40 Ry cutoff energy. The 
equations of motion are numerically integrated within a 
predictor-corrector scheme with a time step of 1.2 fs. 

Since the barrier to associative desorption of hydro- 
gen from Si (100) is rather high (E d « 2.5 eV ||||,||), 
there is no sense in performing molecular dynamics calcu- 
lations starting with the deuterium atoms at the adsorp- 
tion sites because of the extremely low number of des- 
orption events. Therefore we started the desorption tra- 
jectories close to the transition state for dissociative ad- 
sorption which was determined in the earlier study [ff5| . 
The desorbing D2 molecule originates from two D atoms 
which were bonded to the same Si dimer; this pre-pairing 
mechanism is consistent with the observed first-order 
kinetics in experiment [^5|,^6| . In total we have computed 
42 trajectories of D 2 desorbing from Si (100). Eight tra- 
jectories were determined with the Si lattice initially at 
rest, i.e., at a surface temperature of T s = K. Fig- 
ure [l] shows snapshots of such a calculated trajectory. 
The dark Si atoms correspond to the relaxation of the Si 
lattice after the desorption event. Approximately 0.1 eV 
of the potential energy at the transition state is trans- 



ferred to vibrations of the Si lattice which is a rather 
large amount compared to what is known from other 
systems p7| . At the transition state the D-D distance is 
about 40% larger than the D2 the gas-phase bond length; 
consequently molecular vibrations are excited during the 
desorption, as is well known for a long time in associative 
desorption studies (see, e.g., ref. ^Tj). Due to the strong 
anisotropy of the PES, molecular rotations are very ef- 
fectively quenched during the desorption event. 

The mean kinetic energy of D2 desorbing from Si (100) 
was determined by laser induced thermal desorption at a 
rather high surface temperature of T s w 920 K. In order 
to simulate these experimental conditions, we have per- 
formed ab initio molecular dynamics calculations with 
initial conditions corresponding to the experimental sur- 
face temperature. The system was allowed to equilibrate 
for more than 500 fs, whereby the deuterium atoms were 
kept close to the transition state by auxiliary forces F? 
acting on the deuterium atoms 
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The force constant fc = 5.2 x 10 2 Nm" 1 is taken from 
the vibrations of the free hydrogen molecule. The index 
i, i = 1,2,3, denotes the cartesian coordinates, i?* r are 
the transition state coordinates of the two D atoms (see 
Fig. [j]). The extra- force- free region is given by 
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so that the deuterium atoms in z-direction were kept 
closer to the transition state than in x- and y-direction. 

The additional potential was switched off when no ex- 
tra forces were acting on the D atoms, and the energy dis- 
tribution of the D2 molecules desorbing from the thermal 
surface was monitored. Simultaneously trajectories were 
run with the additional potential switched on so that the 
changes in the trajectories of the Si substrate atoms due 
to the desorption event could be followed. This trajec- 
tory was then also used for the next desorption event. 

It is true that thermal desorption usually does not cor- 
respond to a genuine equilibrium situation since it is com- 
monly studied in vacuum. However, desorption events 
are in general not effected by the presence or absence 
of gas molecules in front of the surfaces. This makes it 
possible to treat thermal desorption as if it corresponds 
to a thermal equilibrium situation p8fl . From that it fol- 
lows that we can assume that the phase space density in 
every point in real-space is given by a Boltzmann distri- 
bution, in particular at the transition state. In addition, 
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(E tot ) 


0.72 ±0.17 eV 


{E v ib) 


0.11 ±0.09 eV 


{Ekin) 


0.58 ±0.13 eV 


(Erot) 


0.03 ± 0.05 eV 



FIG. 2. Example of a thermal desorption trajectory which 
was run for 1060 fs. a) x and z coordinates of the trajectories 
of the two deuterium atoms and the two Si atoms of the dimer 
underneath the transition state (TS). The final positions of 
the atoms are marked by the large dots, b) energy redistri- 
bution as a function of the last 140 fs of the run-time of the 
trajectory. Upper panel, solid line: D2 center-of-mass kinetic 
energy, dashed line: kinetic energy of the Si atom closest to 
the transition state (TS); lower panel, solid line: D2 vibra- 
tional kinetic energy, dash-dotted line: D2 rotational energy. 



recrossings of trajectories at the transition state are very 
unlikely to occur in the system D2/Si (100). Due to the 
high binding energy there are large dissipative effects for 
adsorbing molecules; and for desorbing molecules there 
are no restoring forces in the gas phase. This allows 
the application of transition state theory pjJ] and equiv- 
alently also the assumption of thermal equilibrium for 
the initial desorption conditions at the transition state. 

In total 34 "thermal" desorption trajectories were cal- 
culated. Figure |^ shows an example of such a trajectory 
with a total run-time of 1060 fs. For the first 1020 fs the 
molecular dynamics was run with the additional potential 
V c . At t = 1020 fs the extra potential was switched off 
and the molecule was allowed to desorb. The projection 
of the trajectories of the desorbing deuterium molecule 
and of the Si dimer closest to the transition state onto 
the ccz-plane is shown in Fig. ||a). Clearly the vibrational 
excitation of the desorbing D2 molecule can be seen. In 
Fig. |]b) the energy redistribution during the last 140 fs 



TABLE I. Mean energy distribution averaged over 34 tra- 
jectories of D2 molecules desorbing from a Si(100) surface at 
a surface temperature of T a = 920 K (k B T 3 = 0.079 eV). 



of the run is plotted. Besides the vibrational excitation 
the quenching of the rotational motion can be followed. 
Also the acceleration of the D 2 center-of-mass is obvious. 

The mean total, kinetic, vibrational, and rotational 
energies of the D 2 molecules averaged over the 34 ther- 
mal desorption events are listed in table [l| (note that 
k B T s = 0.079 eV at T s = 920 K.). The results show 
vibrational heating, i.e. (£Wf>) > ksT s , and rotational 
cooling, i.e. (E ro t) < ksT s , in agreement with the ex- 
periment [Q. The D2 center-of-mass kinetic energy, 
however, is much larger than the experimental value 
of (E km ) ex P = 0.165 eV > 2k B T s §. The difference 
between the experimental and theoretical results corre- 
sponds roughly to the barrier height = 0.4 eV. A 
closer analysis of the trajectories reveals that still approx- 
imately 0.1 eV of the potential energy at the transition 
state is transferred to the Si lattice, but due to the Si lat- 
tice vibrations the mean adsorption barrier is increased 
by roughly the same amount. Possible contributions to 
the discrepancy between theory and experiment could be: 

(i) insufficient statistics - Certainly 34 thermal trajec- 
tories are not enough in order to claim that a true ther- 
mal average has been performed. However, the widths of 
our distributions seem to be converged quite well and we 
are convinced that by sampling more points of the phase 
space at the transition state as initial conditions the ex- 
cess translational energy of 0.4 eV would not disappear. 

(ii) quantum mechanical effects - There are two im- 
portant quantum mechanical effects not taken into ac- 
count by classical molecular dynamics: tunneling and 
zero-point effects. Low-dimensional quantum dynamical 
studies of D2/Si(100) have already shown that tunneling 
is strongly suppressed due to the thickness of the barrier 
|l^,^8|. As for zero-point effects, it has been shown 
that the sum of all zero-point energies of the hydrogen 
molecule at the transition state approximately equals the 
zero-point energy of hydrogen in the gas-phase. In such a 
situation zero-point effects cancel out as far as the center- 
of-mass kinetic energy is concerned since they lead to an 
almost constant energetic shift of the minimum energy 
paths in the hydrogen coordinates pc[ . Hence we in- 
fer that quantum mechanical effects do not significantly 
change the kinetic energy of desorbing molecules. 

(Hi) dissipation channels not considered - A further 
channel for energy dissipation at surfaces is the excita- 
tion of electron- hole pairs. This channel is not taken into 
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account in the ab initio molecular dynamics simulations 
where the electrons are assumed to follow the ionic mo- 
tion adiabatically in the electronic ground state. While 
electronic excitations of the H2 molecule require a too 
large energy, the excitation of a surface exciton with a 
hole in the dangling bond at the upper Si atom and an 
electron in the lower Si dangling bond is in fact a possi- 
blity. This exciton might have some resemblence to that 
at the (110) surface of GasAs Q The described surface 
exciton at Si (100) will strongly couple to the lattice and 
bring the Si dimer in a geometry parallel to the surface. 
Thus, the H2 desorption ends with an electronic state 
and geometry 0.2 eV higher than the ground state, and 
this energy would be taken out of the H2 kinetic energy. 

(iv) limitations of the GGA functional - It was recently 
shown that the GGA-PW91 functional used in our study 
underestimates the H2 elimination barriers from silanes, 
i.e., small molecules J32|. On the other hand, energies for 
surface systems calculated from cluster models are some- 
times seriously in error (see, e.g., Ref. J33|) and typically 
converge very slowly with cluster size. Hence results for 
small molecules are not directly transferable to surface 
systems. Still the validity of the GGA-PW91 functional 
for the barrier heights of H2/Si(100) certainly remains an 
open question since the application of GGA functionals 
for dissociation at surfaces is still a rather new field. 

(v) experimental uncertainties - The desorption exper- 
iments were done with a Si sample with a high dopant 
concentration of n k 10 19 cm~ 3 . Therefore the Si (100) 
surface used in the experiment might have been metallic. 
It is well possible that for such a system the geometry of 
the clean surface is different from that of an undoped sub- 
strate. Thus, the final state of the desorption is changed 
and possibly also the geometry and energy of the transi- 
tion state. In fact, also in Ref. || it had been speculated 
that high doping might influence the barrier heights in or- 
der to explain that the sticking probability for hydrogen 
molecules at Si (100) is much higher for a highly doped 
sample H than for a weakly doped system pj. 

In addition, recent ab initio calculations indicate that 
laser melting of bulk silicon leads to a transition to a 
metallic liquid state [B4[ which could also occur dur- 
ing the laser-induced desorption used in the experi- 
ment. This would also change the geometrical and elec- 
tronic structure at the transition state. Furthermore, at 
a metallic surface the excitation of electron-hole pairs 
might affect the desorption dynamics and, for example, 
enhance the probability of the above noted exciton state. 

In conclusion, we have presented an ab initio molecu- 
lar dynamics study of the desorption of D2 from Si (100). 
The calculations give a qualified description of the des- 
orption of hydrogen from undoped and good quality 
Si (100). While the energy distribution of the desorb- 
ing molecules in the vibrational and rotational degrees of 
freedom is in agreement with experiment pfl] , the excess 
translational energy is too large compared to the exper- 



iment of Ref. H| . This discrepancy might be caused by 
uncertainties in the determination of the interaction po- 
tential; however, it could also be due to the high doping 
of the Si sample used in the experiment of Ref. || . 
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